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Abstract. The momentum- and frequency-dependent one-body correlation 
function of the one-dimensional interacting Bose gas (Lieb-Liniger model) in the 
repulsive regime is studied using the Algebraic Bcthe Ansatz and numerics. We 
first provide a determinant representation for the field form factor which is well- 
adapted to numerical evaluation. The correlation function is then reconstructed 
to high accuracy for systems with finite but large numbers of particles, for a wide 
range of values of the interaction parameter. Our results are extensively discussed, 
in particular their specialization to the static case. 



One-particle dynamical correlations in the one- dimensional Bose gas 



2 



1. Introduction 

For many decades now, there has been a sustained interest in the physics of one- 
dimensional quantum systems, in no small part due to Bethe's fundamental paper on 
the Heisenberg spin chain [T] . This publication laid the foundation for what has today 
become one of the most important fields of mathematical physics, namely the theory 
of Bethe Ansatz integrable models. 

More recently, renewed motivation for the study of interacting quantum models 
in ID has come from their growing number of experimental realizations. The simplest 
system to consider is probably the one-dimensional gas of bosons with local (delta 
function) interaction (known as the Lieb-Liniger model [2]), which can effectively 
describe atoms confined using optical lattices [3 IH [5l [6l [71 [H HI [TOl [H] . A wide range 
of interactions are accessible, even the strongly-interacting (free fermion-like) limit 
considered long ago by Tonks and Girardeau [12l HI] . 

The integrability of the Lieb-Liniger model is by now textbook material [HI [151 
[TBI [TT] . and many of the standard methods now used in the field were in fact first 
developed and applied to it. For example, after the construction of the ground state 
and elementary excitations in [2j, the Yang- Yang approach to thermodynamics was 
introduced in [18j, leading the way to nonperturbative results on the equilibrium 
properties of the infinite system at arbitrary temperatures. 

The motivation of the present paper is to help further our understanding of 
the interacting Bose gas in one dimension, by addressing one of the long-standing 
unresolved problems associated to it: the calculation of correlation functions. The 
importance of this is clear, in view of the fact that most experimentally-accessible 
response functions are expressible in terms of a few basic dynamical correlation 
functions of local physical fields. 

Because of the relative complexity of the Bethe Ansatz eigenf unctions, it has 
proved much more difficult to provide information on correlation functions than on 
thermodynamics. For the interacting Bose gas, the first set of results probably dates 
back to the work of Girardeau [13j , who considered the strongly-interacting limit. This 
represents a major simplification in that the spectrum becomes that of free fermions. 
The bosonic density operator then coincides with a free fermionic one, in view of its 
blindness to particle statistics. Handling particle statistics in the same limit allows 
to obtain the one-body density matrix as the determinant of a Fredholm integral 
operator, both at zero and nonzero temperatures \T9 \ \2U \ \2T \ 122 ) . 

Away from the Tonks- Girardeau (TG) limit, static local moments of the density 
operator can be obtained. The second moment g2 — {'■ n?{x) :) is given directly 
from the derivative of the partition function with respect to the interaction parameter 
(the Hcllmann-Feynman theorem) [23 , but higher moments g„ with n > 2 are more 
challenging. These were computed for arbitrary temperature in a large coupling 
expansion in [241 125 1 At zero temperature, the third moment 53 was obtained 
exactly for any interaction in [27j using an integrable lattice regularization of the 
Lieb-Liniger model coupled with Conformal Field Theory [551 [Mj considerations. 

The full calculation of correlation functions beyond either the static of local cases 
remains analytically intractable to this day. Besides formulations for the asymptotics 
based on Conformal Field Theory, there exist important results providing determinant 
representations of the exact correlators in terms of auxiliary quantum fields, which 
ultimately map the infinite-volume problem to the calculation of the determinant 
of a Fredholm operator obeying a completely integrable integro-differential equation 
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[5U1 [5T1 [5^ [551 [51] . The asymptotics of the correlation functions can then be 
obtained by solving a Riemann-Hilbert problem (see [TS] and references therein for an 
introduction to the subject). Closed- form analytical expressions for general correlation 
functions have however not yet been obtained. 

Here, we make use of another strategy which one of us originally developed for 
dynamical spin-spin correlation functions in Heisenberg spin chains |351 I36j . Our 
objective is not to give final answers to the theoretical computation of correlation 
functions, but rather to provide a reliable and practical computational scheme for 
obtaining them to very high precision. This will be of great help for experimental 
comparisons (where finite resolutions are inevitable) or comparison with other 
theoretical methods. We specialize here to one-body correlations, following on our 
earlier work |37j on dynamical density-density correlations. 

The overall strategy is summarized as follows: two-point correlation functions 
are expressed as an appropriate sum of form factors over intermediate states, and 
the summation is carried out numerically. The form factors of local operators are 
expressed as determinants of matrices whose entries are simple functions of the sets of 
Bethe rapidities of the two eigenstates involved (for the density, they were obtained in 
[38l [3T] ; see however the more practical representaion we offer here for the field form 
factor). The summation over intermediate states can itself be aggressively optimized 
by making use of the fact that not all intermediate states contribute equally to 
the desired correlator. Such an optimized method is in principle applicable to any 
integrable model as long as the determinant representations of the form factors are 
known, and was baptized the ABACUS method [39j. We here carry on with our 
study of the Lieb-Liniger model by providing the corresponding computations for its 
one-particle dynamical correlation function. 

The plan of the paper is as follows. We first review the basic formulation 
of the theory, and define the correlations we will be interested in. The Algebraic 
Bethe Ansatz for the model is thereafter presented, and used to derive a simplified 
representation for the one-particle form factor in terms of the determinant of a single 
matrix. We then present the results for the correlation functions for a range of 
values of the interaction parameter, and discuss the results in detail, in particular 
their specialization to the static case (momentum distribution function and one-body 
density matrix). Conclusions and directions for future work are presented at the end. 

2. Definitions 

We consider a ring of length L, with N bosonic particles interacting repulsively with 
each other through a local potential of strength c. The time evolution of this system 
is controlled by the Lieb-Liniger Hamiltonian [2] 



An alternate second-quantized formulation is obtained by introducing canonical Bose 
fields ^'(a;) and ^'^(a;) with equal-time commutation relations 




(1) 



[^{x),^Hy)] - Six - y), [^{x),^iy)] = [f^x), ^^y)] = 0. (2) 
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and thus the model corresponds to the quantum nonhnear Schrodinger equation, 

H^f da; [9^*^(x)a^^'(a;) +c^'t(2;)^'t(a;)*(a;)*(a;)] . (4) 
Jo 

The Bethe Ansatz provides a basis for the Fock space. More precisely, in the repulsive 
case c > which we are considering, each eigenfunction in the iV-particle sector is fully 
described by a set of N real parameters {A} at, solution to a set of coupled nonhnear 
equations (the Bethe equations, which we will review and discuss in more details in 
the next section). 

Our aim here is to present computations of the zero temperature one-particle 
dynamical correlation function 

G2(x,<) = (*t(a;,t)*(0,0))jv (5) 

as a function of the interaction parameter c (we consider here only the repulsive case 
c > 0). By identifying the state |{A}) with the ground state and using the closure 
relation of Bethe eigenstates, this correlation function can be represented as a sum 
over form factors of the local field operator between ground- and excited state as 

The evaluation of this expression can be performed in three steps: first, finding 
the solution of the Bethe equations; second, computing the form factors; finally, 
performing the summation over intermediate states. We begin by providing a new 
representation for the form factors, using the Algebraic Bethe Ansatz, which allows 
all these steps to be efficiently implemented numerically for finite numbers of particles. 



3. Algebraic Bethe Ansatz and Form Factors 



The central object of the Algebraic Bethe Ansatz |40] (or Quantum Inverse Scattering 
Method; for an introduction to the concepts and terminology, see [T^) is the R- 
matrix, which solves the Yang-Baxter equation and takes the following form for the 
Lieb-Liniger model (quantum nonlinear Schrodinger equation): 

/ /(m,a) 



in which 



RiX,fi) = 



3(A,m) 





.9(m, A) 
1 




ic 



/(A,A^) 




1 

5(m, A) 


A — /i + i 

X — fi 



\ 



/(m,a) J 



A — /i 

Other standard functions which we will use later are defined as 

X — H -\- ic ^ (ic)^ 



/i(A,^) = 



t(A,/i) 



gi\iJ-) 

h(X,fi) 



(7) 



(8) 



• (9) 



{X — fj.){X — fi -\- ic) 

The monodromy matrix, which will be used to construct all the conserved charges (in 
particular, the Hamiltonian) is represented in auxiliary space as 

■ ^(A) B{X) ■ 
C(A) D{X) 



T(A) 



(10) 
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where the A, B, C, D operators act in the Fock space as foUows. The vacuum is an 
eigenvector of A and D, 

A(A)|0)=a(A)|0), D{\)\0)^d{Xm (11) 

with a(A) = e~*^^/^ and d{\) = e*^'^/^, whereas B and C respectively act as raising 
and lowering operators with the annihilation properties 

(0|B(A)-0, C(A)|0)=0. (12) 

The commutation relations between these operators are quadratic, and given by 

i?(A, ^i){TiX) ® Ti^i)) = {T{y) ® T{\))R{\ /i). (13) 

The Hamiltonian is related by trace identities to the trace of the monodromy 
matrix, A{\) + D{X), and these can therefore be diagonalized simultaneously. The 
(unnormalized) eigenvectors of the Hamiltonian are constructed as 

N N 

\{X}^)=l[B{\,)\0), {{X}n\^{0\1[C{X,), (14) 
i=i i=i 
provided that the rapidities are solution to the Bethe equations 

^.a.l^-qA^^A^^ J = 1,...,A^ (15) 

-^-^ A," — A; — IC 



1^3 

or in logarithmic form, 

1 ^ A - A 27r 
Aj + - ^ 2 arctan^^ - = — 7^, i = l,...,N (16) 

1=1 ^ 

where the quantum numbers Ij are half-odd integers for N even, and integers for 
N odd. Proper eigenfunctions are obtained for sets of non-coincident rapidities 
Aj 7^ Ai, j ^ I. Since the left-hand side of (fTB]) is a monotonous function of Xj, it 
can be proven that all solutions are real [18|; we can moreover span the whole Fock 
space by choosing sets of ordered quantum numbers Ij > > I meaning that 

Xj > Xi,j > I. 

Once the Bethe equations are solved for the rapidities, the eigenstate is fully 
characterized. The state norm is obtained from the Gaudin-Korepin formula [HI HI], 
which can be proven by making use of the commutation relations between B and C 
operators: 

({AV|{AV)=c^ n ^^^detr,g{{X})^\\{X}M\' (17) 

j>k=l 3^ 

in which Q is the Gaudin matrix, whose entries are simple analytical functions of the 
rapidities, 

JV 

L + ^if(A„AO 

1=1 



^jfe({A}Af) = 5jk 
where the kernel is 



-K{Xj,Xk) (18) 



^^^'^)- (A-')'+c2 - ^''^ 
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The energy and momentum of an eigenstate are also given by simple functions of the 
rapidities, namely 

EmN) = -h)^E^, P{{X}n) - E (20) 

j j 
where ft, is a chemical potential used to set the filling in a grand-canonical ensemble. 

For our purposes, we also need explicit expressions for the form factors of the field 
operators '^{x,t) and 'i>''{x,t) between Bethe eigenstates. By making use of the fact 
that the energy and momentum are diagonal operators in this basis, we can explicitly 
take the space and time dependence into account and write 

mN-i\^{x, t)\{\}N) = e^(^--^^)*-*(^--^^)-F({M}^_i; 

({A}jv|*^(:r,t)|{/i}jv-i) =e'(^^-^-^*-'(^^-^-)"i^(MiV-i;{A}A,). (21) 
where 

N-l N 

i^(M^_i;{A}^) = (0| n C{fi,)^{0,0)l[B{X,)\0) (22) 
and F is the conjugate. In these notations, the correlation function ^ becomes 

r,; ll{A}jvirii{Ai}7v--ilr 

in which 

dix, t; A, A*) = (Ex ~ E^)t - i{Px ~ Pp)a;. (24) 

The explicit evaluation of the summation in (j23p to obtain a closed-form 
expression for the correlation function remains an open problem in the theory of 
integrable models. We will here however make use of the building blocks which were 
developed in this spirit. Namely, in [31] . the following useful expression was obtained 
for the field operator form factor: 

F{{^l\N-l■A\^N) = -^^^c 9{l^],lJ'k) Yl 9{>^k,Xj}x 

N-l>j>k>l N>j>k>l 
N N N-l N 

xnnMA;,A,) n d{^,,)l[d{Xl)M^{{X}) (25) 

1=1 j=l J=l 1=1 

Mi{{X}) - ( 1 + 7^ ) detN-iiSjk - aSNk)\a.=o (26) 



d 



bjk - t{fik,Xj) — — - HAj, Mfc)--^— — — ■ (27) 

na=l^(Aa,Aj) na=l'i(Aj,Aa) 



Although equation ((25)) offers an exact representation of the form factors for fixed N, 
(|26p is not convenient for numerical purposes. We therefore here first rewrite the form 
factors in terms of a single matrix determinant without derivative, which is much more 
economically evaluated numerically. The trick is simple, and consists in getting rid of 
the auxiliary parameter a in by using the fact that, given two matrices A and B 
of which B has rank 1, the following identity holds: 

det(A + aB) = (1 - a)dctA + adct(^ + B). (28) 
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Therefore, 

Mi{{X}) = + ^) detN-i{Sjk - aSN.k)\a=o = detN^i{Sjk - ^jvfe). (29) 
Consider now the — 1 x — 1 matrix W with entries 

1 Ua=lifJ-3 - >^a) 



W,^--. -^^W^ —■ (30) 

Since the hnes of this matrix are proportional to -, — -K-^, the determinant of this 



matrix is easily computed: 



N-1 . , N-1 

detw= n n - ^A^)- (31) 

a>htl - a=l 



Consider now rewriting ([^ as 



dct{j:Z-\S,i-Sm)Wik 



detjv-i(5'jfc - SNk) = detW ' ^^"^"^ 

The matrix product above can be calculated explicitly. Consider the sum 

^^c= E^(^''^j)^"^ (33) 
1=1 

and the auxiliary contour integral 

1 / dz{tcr nll(^-Aa) 



I = — (t ^ ' iia=iy (34^ 

27n Ji,\=R^oo (z - \j){z ~ \k)iz - + ic) - ^la) 

The integration contour is a big circle of radius R containing all poles of the integrand. 
The value of this integral is then given by the residue at infinity, which equals zero 
since the integrand behaves like at z — + cxd. On the other hand, the sum of the 
residues at points z — ^a, a. = 1,...,N — 1 inside the contour gives exactly G^^,. In 
addition, there is a pole at z ~ Xj — ic and, for j = k, we also have a pole at z — Xj. 
Since the sum over all of these residues is zero, we obtain 

r+ - ,>A n^#j(^j ~^°) I nf=i(^a-^j-+»c) 

I[a=l (^i - Ma) Xk-X^+ IC J]^^^ (^,^ _ X- + jc) 

Similarly, we have the identity 

N-l 



1=1 



_ Ua^ji^ ~ Ag) ^C X\a=x(Xj »C) 

\\a=l - l^a) ^3-^k+ IC ^^^^ (Aj - + «c). 

This whole procedure allows us to rewrite the matrix product ^Yl!i=\k^fl ^ S]^i)Wik-, 
finally leading to 

M.i[XX))-i -rrjv-i,, ^ detAr_i(7j7c, (^7) 
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in which the matrix C/'s entries are functions of the rapidities of the left and right 
eigenstates, 

UMM, {A}) - ' . ) + ^A^^lStl_ll(^Ki\, - A..) - KiX^ - A,)) (38) 

* lla#j(^a-Aj) 

in which 

na=i(-^a-Aj±zc) 

For the repulsive Bose gas which we are considering, the U matrix has purely real 
entries. The matrix elements also conveniently coincide with those obtained in [38] 
for the (integral of the) local density operator, and which we used in |37| to compute 
density correlation functions. 

4. Dynamical one-particle correlation function 

From the results of the previous section, we can therefore rewrite the one-body 
dynamical correlation function ([5]) as 

G2{x,t)^{^\x,t)^{0,0))N= e^'^"'*^^'''^G(M,{A}), (40) 

{mIn-i 

in which the correlation weight for a given intermediate state is explicitly given by the 
following function of its rapidities {ij-}n-i and of the set of ground state rapidities 
{A}^: 

rdul r^.^ „2A^^l U^.^A^h + [det^-i[/(M, {A})]^ 

where the U matrix entries are given by equation (|38p . and the state norms are given 
by (HID and ([TH]). 

In the present context, it is more convenient to consider the space and time Fourier 
transform of the correlation function, 

G2{k,uj)= dx dte"^'-''''=G2{x,t) 

Jo J -oo 

= 27rL '5(t^-i?w+£^{A})4,P{,,-P{„G({/^},{A}), (42) 

which explicitly displays that each intermediate Bethe eigenstate contributes all its 
correlation weight G at a single point in the fc, lj plane. 

A particular intermediate state can be composed of a number of single-particle 
excitations. Due to the fermion-like structure of the ground state, these take two 
different basic forms, namely (following the terminology of ^) Type I (particle) and 
Type II (hole). Particle-like excitations are obtained by adding an extra occupied 
quantum number outside the ground state interval, while the hole-like excitations are 
obtained by removing one of the occupied ground state quantum numbers. Type I 
modes are Bogoliubov-like particles existing at any momentum, and whose dispersion 
relation (in the infinite system) interpolates between e/(fc) = fc^ at the noninteracting 
7 = point to ei{k) = fc^ -t- 27rn|fc| (where n = N/L is the density) as 7 ^ 00. Type 
II particles do not appear in Bogoliubov theory. Their dispersion relation coincides 
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with the lower threshold of the correlation function, and vanishes at fc = 2fc^, with 
the Fermi momentum given hy kp — nn. Both types of excitations have a dispersion 
relation which approaches fc ^ with a slope equal to the velocity of sound Vs{"f). 

The action of the field operator on a Bethe eigenstate is found to be relatively 
innocuous, in the sense that the size distribution of its matrix elements in this basis 
is very strongly peaked. There are in other words only relatively few form factors of 
large value, and all of these involve only up to a handful of elementary excitations. 
The trace over the Fock space of intermediate states can thus be performed in a highly 
optimized way by searching for and concentrating the computational effort on these 
important intermediate states (part of the idea of the ABACUS method). This is 
implemented with a recursive search algorithm. 

In order to verify the accuracy of our results, we make use of the observation that 
integrating the one-body correlation function (I42p over frequency and summing over 
momenta simply yields back the density of the gas, 

/^SzEG2(fc,-)= g(m,{a}) 

°° fe {m}jv-i 

= G2ix = 0, t = 0) = (^^(0, 0)^'(0, 0)) = n, (43) 

which provides a convenient sumrule for our method. Physically, for large enough N 
and L, there is only one important physical parameter, which is defined as 7 = c/n. 
All the results presented here are obtained by setting the density to one, and varying 
the interaction parameter c. 

4.I. Results 

Our results for the dynamical one-body correlation function ([5]) for different values 
of 7 are presented in Figures [T] and [2] as density plots. To obtain continuous curves, 
the energy delta functions in ((42|) are smoothencd into Gaussians of width slightly 
larger than the typical interlevel spacing. The color coding of these figures follows 
the logarithm of the correlation function. All data presented in these was obtained by 
considering systems of = 150 particles with unit density. For each individual plot, 
around seventy million intermediate states were taken into account, yielding extremely 
good saturation of the sum rule (|43|) , as summarized in Table [1] The summation we 
actually perform is restricted to intermediate states within a window of momentum 
fc G [0,4fcF], composed of (possibly many) single-particle excitations of momentum 
1^1 < kmax = Skp- The correlation weight carried by states outside of this sub- Fock 
space is negligible, as can be clearly seen from the sum rule saturations achieved. This 



Table 1. Sum rule saturation percentages achieved as a function of the interaction 
parameter, for various values of TV. The last line gives the contribution coming 
from states outside of the sub-Fock space used for the numerics, as computed 
from the asymptotic form of the static correlator (see text). 
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Figure 1. Pseudocolor plots of the logarithm of the dynamical one-body 
correlation function of the Lieb-Liniger gas, for 7 = 1/4,1/2,1,2,4 and 8. The 
data was obtained from systems at unit density and A'^ = 150 particles. The 
horizontal axis is momentum, and the vertical one energy. Sum rule saturations 
are listed in Table [T] At small 7, the correlation weight at fixed momentum falls 
mostly around the Type I dispersion relation. 
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Figure 2. Same as Figure [T] but now for 7 = 16,32,64 and 128, showing the 
progression of the distribution of correlation weight in the strongly-interacting 
regime. As 7 increases, more correlation weight gets concentrated near the lower 
threshold. 



can also be independently checked because the large k tail of the static correlation 
function is exactly known from [42j (as discussed more extensively in the next section, 
see Eq. (j49l) below). In the last row of Table [T] we report the contributions calculated 
from Eq. for all the states with k > Akp. We checked that the error due to 

finite and to subleading corrections to Eq. (|49)) do not change the accuracy we 
reported. From these data it is evident that states with k > 4kp contribute for almost 
all the missing part for = 50 and for a relevant part for N = 100. In all cases, 
the (rest of the) missing part of the sum rules is therefore ascribed to intermediate 
states within the sub-Fock space we consider, but which we do not include in the 
summation. In fact, our sum only considers an extremely small fraction of all the 
available intermediate states within the sub-Fock space, made of only a handful of 
elementary excitations. The distribution of correlation weight between 0,2,4,6 and 
8 particle states for a given value 7 = 16 for different values of N is given in Table 
2 (here, 2n particle states are made from n particle-hole excitations on the ground 
state of — 1 particles). These numbers represent lower bounds of the contributions 
from the different classes of states. There are large variations of this distribution as a 
function of N, although the full correlator itself remains more or less invariant. The 
general scenario is that and 2 particle contributions are strictly decreasing functions 
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Table 2. Distribution of sum rule weight between intermediate states involving 
Np particles, with Np = 0, 2, 4, 6 and 8. The numbers given are the fraction of 
the total sumrule associated to each sub-type of excited state which was obtained 
during the runs. We present data for 7 = 16, for three values of A'^. These 
are lower bounds only, since we only scan over a very small fraction of the total 
number of multiparticle states. 
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of N, while higher particle number contributions first increase and thereafter decrease, 
with a maximum at N^''^ such that > if n > m. As a function of 7, 

we observe as naturally expected that for small 7, intermediate states with smaller 
(larger) numbers of particles carry more (less) weight, with more weight shifting to 
higher particle numbers when 7 increases. 

We also provide in the left panels of Figures 3 and 4 a series of fixed-momentum 
cuts for a number of values of the interaction parameter, with data for two different 
system sizes, smoothened as explained above (without smoothing, these would be 
6 peaks; the smoothing however blurs the lower threshold slightly). At small 7, 
the correlation is peaked around the Type I dispersion relation. As the interaction 
parameter increases, the peak becomes broader and moves towards lower energies, as 
expected in view of the increasing fermion-like nature of the system. At intermediate 
values of fc, the correlation width also increases with 7, with the peaks eventually 
shifting altogether from high to low energies. The difference between the N — 100 
and iV = 150 curves can mostly be ascribed to imperfect sumrule saturation (the 
logarithm of the correlation function is plotted, and therefore the actual difference 
between curves lies in regions where the correlation weight is extremely small) . 

On the right panels of these figures, the same data sets are used to plot the 
integrated correlation, 

Gi{k,uj) = / duj'G2{k,uj) = 
Jo 

= 27rL e(c^'-i?|^|-|-£;{A})4,p,,,-p^.jG({/i},{A}), (44) 

for which no smoothing procedure is necessary. These curves are in fact series of steps 
of height corresponding to the correlation weight of any intermediate state lying at 
this point in fc, ui space. In the infinite size limit, these would become smooth curves 
whose oj derivative would be the one-body correlation function. The discreteness of 
the steps due to the finite number of particles can just be made out in the plots. The 
fact that no artificial smoothing of the delta functions is used in these data sets makes 
them the most appropriate objects for eventual comparison with other theoretical or 
numerical methods (as e.g. those employed in Refs. [43! for the dynamical structure 
factor). 
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a/kt 



Figure 3. Left; fixed-momentum cuts of the logarithm of the one-body 
correlation function, for 7 = 1/4, 1 and 4. As 7 increases, the peaks widen. 
Right: logarithm of the integrated correlation function (144 l l for the same values 
of k and 7 (e = 10~^^ is a regulator for the logarithm). 
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Figure 4. Same as above, now for 7 = 16 and 64. The peaks have moved to low 
energy. 



5. Momentum distribution function and one-body density matrix 

From the data of the dynamical one-particle correlation function presented in the 
previous section it is very easy to recover the static correlation function 

/>00 7 

n(fc)- ^G2ik,Lu). (45) 

n{k) represents the relative occupation of the single particle state of momentum k and 
is consequently known as the momentum distribution function. It is the quantity that 
is directly measured in the ballistic expansion of a trapped gas (see e.g., the review 
|44j ) . Its Fourier transform is the one-body density matrix 

p{x)^ 'M^ ron~^ 



L \ L 

J^^ + X,X2,-.-, XN)'^oixi,X2, . . . , XN)dX2 ■ ■ ■ dxN 

Jq |^'o(a:;i,a;2, . . . ,XNWdxi . . . dxN 



,(47) 

where "i'o{xi,X2, ■ ■ ■ ^x^^) is the A^-body ground-state wavefunction. p{x) is the 
reduced density matrix of a single particle when all the other degrees of freedom 
have been integrated out. 
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Figure 5. Momentum distribution function for some representative value of 
7. Left: n(k) in log-log scale. In the inset the large momentum section of and 
the asymptotic law. Right: plot of kn(k) as function of k in log-scale (this 
can be conveniently compared with the data of Ref. [55]). The exactly known 
Tonks-Girardeau (7 — > 00) limit is plotted as reference. 



Given its theoretical and experimental importance, the momentum distribution 
function is the correlation function that has been mostly studied in the literature. In 
the TG limit, its exact expression is known since the work of Lenard [19j : for finite N 
and L, p{x) can be written as a simple Toeplitz determinant (see also |45j ) and the 
momentum distribution function can be readily obtained via Fourier transform. The 
Lenard formula has been re-obtained in the framework of algebraic Bethe Ansatz in 
Ref. [22j . Some analytical asymptotic expansions in the TG limit have been reported 
in Refs. [46 l [2T | l45 t l47] and a I/7 expansion has been developed in [48]. We mention 
that these results have been recently generalized to a lattice of impenetrable bosons 
[49] and to anyonic gases [50| . 

In the case of finite coupling constant, only the asymptotic expansion for small 
and large momenta are known. The infrared behavior can be obtained from the 
hydrodynamic theory of low-energy excitations [51j 

n(k) cx k°''^ , (48) 

where a = Vs/ (47rn) and Vg is the speed of sound. The above result holds for k <C S,'^ , 
where ^ = ^/2/vs is the healing length. The corresponding behavior of the one-body 
density matrix is p{x) oc x~°' for x ^ ^. The exponent a interpolates between and 
1/2 when the coupling constant 7 is increased from to 00. The large momentum tail 
can be elegantly obtained by using simple theorems on the Fourier transform [42| and 
can be related to the previously calculated second moment of the density operator 
[24], yielding {n{k) corresponds to w{p) in Ref. [42] and not to W{p)) 

Nn{k)=-f^e'{-f)(J^y , for fc ^ 00 (49) 

where 6(7) is the ground-state energy per particle in the thermodynamic limit (in the 
TG limit this formula was previously obtained in Ref. [47[ 145] ). Similarly a small x 
expansion can be written for the one-body density matrix: 

P^^^l + yc^\nx\\ (50) 

4=1 
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The fc~'^ law for n{k) implies that the lowest odd power in this Taylor expansion 
is thus Ci = 0. The coefficient C3 is easily read from Eq. yielding 

C3 = 7^e' (7)712 C2 can be obtained [H] using the Hellmann-Feynman theorem: 
C2 = — l/2[e(7) — 7e'(7)]. The other coefficients Ci with i > 3 are still unknown. 

To the best of our knowledge the results reported above are all that is exactly 
known about the momentum distribution function (some results for finite and small 
A'^ are reported in Ref. [52j). In particular nothing is known for finite coupling in 
the intermediate window of momenta, that is actually the region easier to access in 
experiments (in fact the infrared behavior is usually obscured by the harmonic trap, see 
e.g. [7l[53l[54l[55l[56]). As a consequence, up to now this interesting regime has been 
investigated only with purely numerical methods as Quantum Monte Carlo [55j and 
density matrix renormalization group [53]. However, these methods, although very 
powerful, were unable to describe all the physics from the infrared to the ultraviolet. 
Our results demonstrate that the ABACUS method is able to do this. 

In Fig. [5] we present n{k) as function of k for several values of the interacting 
parameter and for N — 100, 150. The curves clearly show that finite size effects are 
negligible on this scale. We also plot kn(k) for comparison with the Quantum Monte 
Carlo data of Ref. [55]: the data shows an overall qualitative agreement. The exact 
result of Lenard for the TG limit is also given for comparison. 

One advantage of our method is that we are able to reproduce the k~'^ law to 
high accuracy. In the inset of Fig. [5] we plot n{k) versus k for large momenta and for 
several values of 7. The k^^ tails are evident, and agree with the asymptotic curves 
predicted from Eq. (|49]) . We stress that in these curves there is no fitting parameter, 
both the power-law and the numerical prefactor are fixed according to Eq. (|49]) . This 
therefore confirms that our method provides accurate results even at high momentum. 

Now we turn to the real spatial dependence of this correlation function. From our 
data, the one-body density matrix is easily obtained by means of Fourier transform 
according to Eq. (|47]) . In Fig. [6] we plot p{x) as function of x/L for some values of 
7. The inset in log-log scale clearly shows the power-law behavior x~°' for a; 3> ^ 
(but smaller than L/2). At large distances there are evident deviations due to 
the finiteness of the system. Conformal field theory [S7] predicts that for periodic 
boundary conditions the correlation function gets modified according to 



For X this form reproduces the power-law but it also completely describes 

the behavior for any x with the only condition x 3> ^. In fact, plotting L°'p{x) 
versus x/L all the curves for various values of N fall on the same master curve, which 
depends on 7 only through the exponent a. The proportionality constant in Eq. (|5ip 
can be fixed by imposing a given value ai x — L/2, but in Fig. [6] we let it free and 
consequently there is no fitting parameter again. It is evident that only for very small 
value of x/L it is possible to notice the effect of the finite healing length ^. 

The small x behavior is also easily accessed with ABACUS. In Fig. [7] we report 
p{x) for very small x (up to a; = 1) for N = 100. Compared to Fig. [5] we are now 
showing a very small fraction of the total plot. The numerical data are compared 
with the Taylor expansion Eq. (|50p with the coefficients C2 and C3 evaluated in Ref. 
[42] . For very small x the agreement is perfect and again there is no fitting parameter. 
However, we notice that Eq. (|50p with only two c^'s describes a very limited region, 
probably very difficult to reach in real systems. A curious fact can be easily deduced 




(51) 
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Figure 6. 

One-body 
density matrix 
p(x) for 7 = 
64,8,1,0.25 
and N = 
100, 150, 200 
compared with 
the finite-size 
scahng (FSS) 
form. Inset: 
The same plot 
in log-log scale 
to show the 
power law 
behavior for 
€ « X < L/2. 



from Fig. [6l the correction to the asymptotic result is positive for large 7 and negative 
for small values, so the unknown coefficient C4 changes sign as 7 increases (contrary to 
C2 and C3, which are always respectively negative and positive). From our data, this 
change should occur around 7^8. 

6. Conclusion 

In conclusion, we have presented results for the one-body dynamical correlation 
function of the one-dimensional Bose gas (Lieb-Liniger model) by using the ABACUS 
method, mixing integrability and numerics. We obtained the full momentum and 
frequency dependence of the correlations, for values of the interaction parameter 
interpolating continuously from the weakly-interacting regime to the strongly- 
interacting Tonks-Girardeau regime. We wish to stress that the present method not 
only yields previously inaccessible results on the dynamics, but by extension also allows 
to go beyond other methods when restricting to static quantities. For example, we 
have showed how our results, integrated to recover the static correlation function, are 




Figure 7. Small x behavior of the one-body density matrix p{x), compared with 
the asymptotic result Eq. I|50| l. Left: 7 = 0.25,0.5,1,2 from top to bottom. 
Right: 7 = 4, 8, 16, 32, 64 from top to bottom. 
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the first to be in agreement with all asymptotic predictions. 

Wc hope that our results will provide motivation for further experimental work 
on quasi-one-dimensional gases, in particular on the direct observation of dynamical 
correlation functions in these systems. On the other hand, we intend to use the 
correlation fimctions we have obtained for the idealized Bose gas as a starting point for 
addressing more general situations, including the effects of intertube coupling and/or 
confining potential. Other possibilities include the study of entanglement in this 
stongly-correlated system, using measiires expressible in terms of simple combinations 
of the correlation functions we have studied. We will report on these issues in future 
publications. 
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